function[IRF]=LP_Jorda_model(data,hor,p,identification,scale)

%----------------- LP (Jorda)------------------
Y=data;
[~,rhs,lhs]=lag_variable(Y,p);
[T,N]=size(lhs);
M=size(rhs,2);

%Estimate LP at horizon 0 and 1 (same as VAR)
Phi=(rhs'*rhs)\(rhs'*lhs);
e=lhs-rhs*Phi;
sigma=e'*e/(T-M);

Phi_int=Phi(2:end,:);
F_lp = [Phi_int';eye(N*(p-1)) zeros(N*(p-1),N)];
bigeye= [eye(N) zeros(N,N*(p-1))];


%% Identification
% identification=1 if Cholesky
% identification=2 if Long-Run

if identification==1 % Cholesky
    if scale==0
      A0=chol(sigma,'lower');
    elseif scale==1
      A0=chol(sigma,'lower');
        for i=1:N
         A0(:,i)=A0(:,i)./A0(i,i);
        end  
    end
elseif identification==2 % Long-run
   invC =  eye(p*N)-F_lp;
      C =  inv(invC);
      C =  C(1:N,1:N);
      D =  chol(C*sigma*C','lower');
     A0 = inv(C)*D;      
end

%% IRF estimation
IRF=zeros(N^2,hor+1);

%----- For h=0,1-----
for j= 1:2
    if j==1
    IRF(:,j)= reshape(A0,[],1);  

    elseif j==2
    IRF(:,j) = reshape(bigeye*F_lp*bigeye'*A0,[],1);
 
    end
end


%----- For h>1-----
for j=3:hor+1
    rhs_h=rhs(1:end-j+2,:);
    lhs_h=lhs(j-1:end,:);  
    Phi_h=(rhs_h'*rhs_h)\(rhs_h'*lhs_h);
    Phi_h_int=Phi_h(2:end,:);

% Companion matrix
Fh_lp = [Phi_h_int';eye(N*(p-1)) zeros(N*(p-1),N)];

% IRFs
IRF(:,j)=reshape(bigeye*Fh_lp*bigeye'*A0,[],1);

end



end
